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FOREWORD 


The  work  reported  herein  was  conducted  by  General  Applied  Science  Labora¬ 
tories,  Inc.,  principally  under  contract  F3361 5-78-C-301 6  issued  by  the 
U.S.  Air  Force,  Air  Force  Systems  Command,  Air  Force  Wright  Aeronautical 
Laboratories,  Air  Force  Flight  Dynamics  Laboratory,  Wright-Patterson  AFB, 
Ohio  45433,  under  the  technical  direction  of  Dr.  George  Seibert,  AFWAL/FIME. 
The  work  was  also  partially  supported  by  GASL  internal  funds. 

The  report  consists  of  two  parts.  The  first  volume  describes  the  theoreti¬ 
cal  formulation  and  method  of  analysis  employed  in  the  subject  computer 
program.  The  second  volume  is  a  guide  to  the  use  of  the  program. 

The  cooperation  of  Dr.  John  Lordi  of  CALSPAN  Corporation  in  furnishing  a 
listing  of  their  one-dimensional,  nonequilibrium  streamtube  program  and 
related  reports  is  gratefully  acknowledged.  Portions  of  the  CALSPAN  code 
have  been  incorporated  as  subroutines  of  the  subject  program,  thereby 
providing  a  degree  of  commonality  with  other  nonequilibrium  flow  programs  in 
current  use  in  the  Thermomechanics  Branch  of  the  Air  Force  Flight  Dynamics 
Laboratory. 


The  work  was  performed  during  the  period  15  May  1978  to  30  October  1979. 
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INTRODUCTION 


The  Fortran  Program  R2D2  is  designed  to  carry  out  a  finite  difference  solution 
of  the  unsteady  equations  of  motion  for  an  inviscid  chemically  reacting  and/or 
Vibrational ly  excited  mixture  of  gases,  in  a  two-dimensional  (planar  or  axisym- 
metric)  duct.  The  program  is  written  in  Fortran  IV  for  use  on  a  CDC  6600 
computer  operating  under  NOS.  Two  versions  of  R2D2  have  been  developed:  one 
for  nonequilibrium  chemistry  and/or  vibrational  relaxation,  R2D2NE,  and  the 
other  for  fully  equilibrated  chemistry  and  vibrational  excitation,  R2D2EQ.  The 
present  version  of  R2D2NE  requires  256K  (octal)  words  of  core  storage,  and 
R2D2EQ  requires  15**K  (octal).  An  additional  56K  (octal)  words  are  required  to 
compile  the  source  codes.  The  FTN  compiler  at  0PT=1  has  been  used  for  the  test 
cases,  although  a  more  efficient  executable  code  may  result  from  use  of 
0PT=2.  However,  in  the  authors  experience,  some  caution  must  be  exercised  with 
due  to  occasional  compiler  errors. 

The  present  report  describes  the  input  and  output  formats  of  the  programs,  the 
operational  procedures,  a  synopsis  of  the  Main  Program,  principal  subroutines 
and  utility  subroutines,  examples  of  the  output,  and  related  information  pertinent 
to  the  use  of  the  program. 
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GENERAL  DESCRIPTION 


1.  INPUT  DATA  REQUIREMENTS 

The  required  information  which  must  normally  be  provided  as  input  data  consists 
of: 

a.  Duct  geometry 

b.  Boundary  conditions  (inlet/exit  stations) 

c.  Initial  conditions 

d.  Finite-difference  grid  specifications 

e.  Output  format  options 

The  duct  geometry  is  defined  by  a  series  of  arbitrarily  spaced  coordinates 

along  the  outer  wall,  (RTIP,  ZTIP),  and  along  the  inner  wall,  axis  or  plane 

of  symmetry  (RHUB,  ZHU,B) .  (Note  that  the  Fortran  names  generally  refer  to 

a  cylindrical  coordinate  system,  i.e.,  z=x^ ,  r^x^ ,  vz=u^»  vr=u2’  ve=u3’  etc 

However,  in  planar  cases,  they  relate  to  a  Cartesian  system,  e.g.,  z=x, 

r=y,  v  =v  ,  v  =v  ,  v  =0,  etc.)  The  inlet  and  exit  planes,  as  well  as  a 
'  zx  ry  B 

series  of  planes  which  define  the  boundaries  of  computational  domains  into 
which  the  duct  may  be  divided,  are  defined  by  pairs  of  axial  coordinates  ZBDH 
and  ZBDT  representing  the  intersection  points  of  each  plane  with  the  duct 
walls.  This  data  is  shown  schematically  in  Figure  (1). 

The  boundary  conditions  consist  of  certain  combinations  of  flow  properties 
at  the  inlet  and  exit  stations,  as  described  in  Volume  I  of  this  report. 

These  may  be  given  at  a  series  of  arbitrarily  spaced  radial  positions  along 
the  inlet  and  exit  planes,  given  by  RIN  and  ROUT,  respectively.  Since  the 
specified  values  at  these  positions  will  be  interpolated  and/or  extrapolated 
to  the  inlet  and  exit  grid  points,  the  first  and  last  values  of  RIN  and  ROUT 
need  not  lie  on  the  duct  walls  or  axis.  See  Figure  (2). 

The  data  to  be  specified  at  the  inlet  consists  of  (a)  either  total  pressure, 
PTINLI,  and  total  temperature,  TTINLI,  or  static  pressure,  PINLEI,  and  static 
temperature.  TINLEI,  (b)  either  radial  velocity,  VRINI,  or  radial  flow  angle, 
ALFA  I ,  or  ratio  of  c? rcumferent ia I  vorticity  component  to  density,  WETAI,  and 
(c)  either  circumferential  velocity,  VTH1NI,  or  circumferential  flow  angle,  BETAI, 
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or  ratios  of  circumferential  and  axial  vorticity  components  to  density,  WXII 
and  UZETAI.  The  inlet  axial  velocity  component  VZINI  is  also  specified, 
however  it  is  only  used  if  the  local  Mach  number  is  supersonic.  At  the 
exit  boundary  either  the  total  pressure,  PTOUTI  or  the  static  pressure, 

POUTLI  may  be  specified.  If  the  exit  axial  velocity  is  supersonic,  the 
pressure  boundary  condition  will  be  automatically  replaced  by  a  solution  of 
second  characteristic  compatibility  relation.  Although  use  of  the  total 
pressure  as  the  exit  boundary  condition  is  admissible,  it  has  been  found  to 
be  unstable  and  is  therefore  not  recommended. 

In  addition,  the  species  concentrations  ALPHAI  are  specified  at  the  inlet 
for  the  nonequilibrium  version,  R2D2NE.  In  the  equilibrium  version,  R2D2EQ, 
the  values  of  ALPHAI  denote  elemental  concentrations. 

The  initial  conditions  are  obtained  in  the  manner  described  in  Volume  I  of 
this  report  using  a  specified  reference  or  plenum  total  pressure,  PTINF, 
total  temperature,  TTINF,  and  Mach  number,  EMINF.  A  corresponding  static 
pressure,  PINF,  temperature,  TINF,  and  sound  speed,  AINF,  are  computed  from 
the  specified  reference  or  plenum  data,  and  used  to  nondimensional i ze  all 
flow  variables  used  within  the  main  body  of  the  program.  A  reference  length 
REF  is  also  used  to  nondimensional ize  all  lengths. 

The  finite  difference  grid  is  defined  by  selecting  the  number  of  grid  col¬ 
umns,  JMAX,  and  grid  rows,  KMAX,  for  each  domain.  The  domains  extend  from 
the  lower  wall  or  axis  to  the  upper  wall,  and  from  (a)  inlet  station  to  the 
boundary  plane  of  the  first  domain  designated  by  ZBDH(2)  and  ZBDT(2) ; 

(b)  from  there  (ZBDH(2) ,  ZBDT(2))  to  the  next  boundary  plane  (ZBDH(3),  ZBDT ( 3) ) ; 

(c)  etc.,  to  the  exit  plane,  which  terminates  the  last  domain.  In  general,  a 
domain  extends  from  the  boundary  plane  through  (ZBDH(l),  ZBDT(l))  to  the 
boundary  plane  through  (ZBDH(I+1),  ZBDT { 1+ 1 ) ) .  The  total  number  of  domains 
is  specified  by  IDOM. 

The  normal  output  format  consists  of  tabulations  of  flow  properties  at  each 
grid  point,  at  selected  time  intervals;  this  output  format  is  arranged  by 
grid  columns.  Provisions  for  output  along  specified  stations  and  along 
specified  streamlines  is  also  included  in  the  program,  but  the  associated 
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coding  for  printout  of  the  thermochemi cal  data  has  not  been  completed. 

The  data  which  defines  the  chemical  components  of  the  mixture  of  gases,  their 
thermodynamic  properties  and,  in  R2D2NE,  their  chemical  reaction  kinetics, 
chemical  rate  constants  and  vibrational  relaxation  rates,  are  contained  in 
a  Block  Data  Subroutine.  In  R2D2EQ  the  chemical  and  vibration  rate  data 
are  replaced  by  equilibrium  composition  data.  The  present  version  of 
R2D2NE  includes  data  for  an  8  species,  10  reaction  model  of  high  temperature 
air,  taken  directly  from  References  (1-3).  Consideration  of  other  chemical 
systems  will  require  recoding  of  the  Block  Data  Subroutine,  and  corresponding 
reidentification  of  the  species  labels  in  the  input  and  output.  If  the 
system  includes  more  than  8  species  or  more  than  10  reactions,  correspond¬ 
ing  increases  in  the  dimensions  of  certain  variables  will  also  be  required. 

In  view  of  the  large  size  of  the  present  version  of  R2D2NE,  the  dimensions 
have  been  kept  at  the  minimum  needed  for  the  air  model. 

It  is  also  pointed  out  that  the  equilibrium  version,  R2D2EQ,  employs  the  same 
species  model  of  air  as  defined  in  R2D2NE.  It  has  been  assumed  that  the 
equilibrium  concentrations  of  the  8  species  are  only  a  function  of  pressure 
and  temperature,  which  implies  that  the  elemental  fractions  are  constant. 

This  would  be  true  for  a  uniform  distribution  of  the  elements  across  the 
inlet  station.  However,  if  the  elements  were  nonuniformly  distributed 
across  the  inlet,  then  the  elemental  fractions  would  differ  from  stream¬ 
line  to  streamline.  In  this  case  the  equilibrium  data  would  have  to  include 
the  elemental  fractions  as  parameters.  In  view  of  this  possibility,  R2D2EQ 
includes  the  solution  of  species  conservations  equations  for  the  elemental 
concentrations.  The  current  version  is  restricted  to  a  binary  system  of 
elements,  which  are,  as  noted  above,  irrelevant  to  the  equilibrium  concen¬ 
trations  of  species  in  the  assumed  model. 

Careful  attention  should  be  paid  to  the  system  of  units  for  the  input  data. 
All  input  (except  the  Block  Data  information)  is  nondimens ional i zed  in 
Subroutine  INDATA,  and  all  output  data  is  consistently  red i mens i ona 1 i zed 
in  Subroutine  OUTPUT.  Thus,  any  self-consistent  set  of  units,  e.g.,  SI, 
conventional  English,  CGS,  etc.,  can  be  used,  subject  only  to  the  restriction 
that  it  must  also  be  consistent  with  the  units  used  for  the  thermochemical 


data  in  the  Block  Data  Subroutine.  The  spatial  coordinates  present  an  ex¬ 
ception,  however,  since  they  need  only  match  the  units  of  the  length  scale 
REF,  with  respect  to  which  they  are  nondimensional ized,  and  REF  may  be  in 
any  set  of  units  (or  may  be  nondimensional),  independently  of  the  other  data. 

The  current  version  of  R202  employs  Block  Data  adopted  from  References  (1-3). 

Therefore,  all  input  data  (except  the  spatial  coordinates)  must  be  in  the 
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CGS  system,  e.g.,  pressure  in  dynes/cm  ,  temperature  in  K,  velocity  in  cm/sec, 
etc. 
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2.  PROGRAM  OPERATION 

Any  particular  run  should  proceed  through  a  specified  number  of  time  steps, 
from  an  initial  time  count  of  I T I  ME  to  a  final  time  count  of  ITIMEF.  Deter- 
miniation  of  whether  an  asymptotic  (steady)  solution  has  been  reached  at 
ITIMEF  must  be  made  by  examining  the  solution  at  several  intermediate  time 
steps  preceding  ITIMEF.  Since  the  number  of  time  steps  required  to  achieve 
a  steady  solution  may  not,  in  general,  be  accurately  predictable,  the  capa¬ 
bility  to  restart  a  run  has  been  provided.  The  restart  option  flag  ISTART 
designates  whether  a  run  is  the  initial  statement  of  the  problem  (ISTART=0) 
or  a  restart  of  a  previous  run  (l$TART=1  or  =2).  A  restart  run  with  ISTART=2 
represents  a  virtually  uninterrupted  continuation  of  a  previous  run,  with 
only  a  redefinition  of  the  initial  and  final  time  counters,  ITIME  and  ITIMEF, 
and  a  possible  alteration  of  the  output  options.  A  restart  run  with 
ISTART=1  can  utilize  the  previous  flow  solution  with  a  redefinition  of  duct 
geometry  and/or  boundary  conditions. 

Magnetic  tape  or  disc  files  designated  TAPE1  and  TAPE2  are  used  as  the  input 
media  for  a  restart  run,  and  TAPE3  and  TAPE4  are  used  as  the  output  media 
for  storage  of  data  at  termination  of  a  run  for  possible  subsequent  use  in 
a  restart  run  (i.e.,  TAPE3  and  TAPEA  are  redesignated  as  TAPE1  and  TAPE2) . 
However,  these  files  are  not  used  during  the  course  of  a  run. 


3.  OUTPUT  DATA 

The  frequency  at  which  output  is  obtained  at  the  grid  points  is  controlled 
by  the  input  parameter  IMESHP,  the  frequency  at  which  streamline  output  is 
obtained  is  controlled  by  ISTREP,  and  the  frequency  of  station  output  is 
controlled  by  ISTATP.  Grid  point  output  can  be  suppressed  by  setting  IMESH=0, 
streamline  output  can  be  suppressed  by  setting  ISTREM*0,  and  station  output 
can  be  suppressed  by  setting  ISTAT=0.  (The  same  effect  can  also  be  obtained 
be  obtained  by  setting  IMESHP,  ISTREP  or  ISTATP  to  a  value  larger  than 
ITIMEF.)  The  output  data  is  summarized  as  follows.  It  should  be  noted  here 
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♦ 


that  although  the  ISTREM  and  ISTAT  options  have  been  Included  in  the  input, 
they  have  not  been  fully  coded  in  the  present  versions  of  R2D2  and  their  use 
is  discouraged. 

The  grid  point  output  consists  of  the  values  of  the  three  components  of 
velocity,  VZ,  VR  and  VTH,  the  static  pressure  and  static  density,  P  and  RHCf, 
the  static  temperature  T,  and  Mach  number  MACH,  the  weight  flow  rate  MOOT 
and  the  entropy  ENTRO.  The  second  block  of  output  are  the  species  array 
ALPHA.  In  the  present  versions  of  the  program.  Block  Data  Statements  define 
the  order  of  the  species  as  follows: 

1st  column  -  oxygen  atom  (0) 

2nd  column  -  nitrogen  atom  (N) 

3rd  column  -  electron  (e") 

4th  column  -  argon  (Ar) 

5th  column  -  oxygen  molecule  (0^) 

6th  column  -  nitrogen  molecule  (N^) 

7th  column  -  nitricoxide  molecule  (NO) 

8th  column  -  nitricoxide  ion  (N0+) 

The  third  block  of  output  for  the  nonequilibrium  version  consists  of  the 
vibrational  temperatures  for  and  N2>  respectively.  In  the  equilibrium 
version,  R2D2EQ,  this  block  also  includes  the  elemental  mass  fractions. 

The  fourth  block  of  output  consists  of  the  mixture  static  enthalpy,  ENTH, 
mixture  specific  heats  (per  unit  mass)  including  all  vibrational  degrees 
freedom1*, CPV  and  CVV,  the  mixture  specific  heats  (per  unit  mass)  excluding 
the  nonequilibrium  vibrational  contributions,  CP  and  CV,  the  ratio  of 
specific  heats  (y^),  GAMF,  and  the  equilibrium  isentroplc  exponent  (yfi) , 
GAME  (in  R2D2EQ  only). 

(It  should  be  noted  that,  by  definition,  the  specific  heats  CP  and  CV  will 
not  include  the  vibrational  contributions  even  as  TVj->T,  since  the  condition 
TypT  is  not  explicitly  recognized  by  the  program.  Thus,  as  TVj-»-T  for  all 
i  species  in  the  range  f < I <jg ,  CPV  and  CVV  acquire  the  physical  significance 
usually  attributed  to  specific  heats  at  constant  pressure  and  volume, 
respectively.) 
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The  system  of  units  in  which  the  variables  are  printed  will  correspond  to  that 
in  which  the  input  is  stated,  except  for  the  last  block  which  has  no  counter¬ 
part  in  the  input  data.  The  static  enthalpy  ENTH  will  be  in  the  units  of 
veloci ty- squared,  and  the  specific  heats  will  be  in  units  corresponding  to 
the  species  specific  heats  given  in  the  Block  Data  Subroutine. 

The  streamline  output  (should  it  be  activated)  consists  of  linearly  inter¬ 
polated  data  from  the  flow  solution  at  the  grid  points.  The  variables  printed 
are  the  magnitude  of  the  velocity  vector,  V,  the  radial  and  circumferential 
flow  angles,  ALPHA  and  BETA,  the  static  pressure,  P,  static  density,  RH(f 
and  the  ratio  of  the  difference  between  the  mixture  entropy  and  the  reference 
or  plenum  entropy  to  the  gas  constant,  ENTRO.  ALPHA  is  defined  as  the  angle 
between  the  radial  and  axial  velocity  components,  and  BETA  is  the  angle  be¬ 
tween  the  circumferential  and  axial  components.  The  entropy  is  defined  using 
input  value  of  y,  and  neglecting  entropy  due  to  nonequilibrium  reactions. 

The  station  output  also  consists  of  linearly  interpolated  data  from  the  flow 
solution  at  the  grid  points.  The  variables  printed  are  the  same  set  identi¬ 
fied  above  in  the  streamline  outputs. 

Sample  printouts  from  both  versions  of  R2D2  are  included  in  the  Appendix. 


LIST  OF  INPUT  DATA  FOR  PROGRAM  R2D2NE  (NONEQUILIBRIUM) 


Card 


No. 

Columns 

Format 

Description 

1 

1-80 

Hollerl th 

Title  card  to  Identify  run. 

2 

1-15 

£15. 0 

TIME-Elapsed  time  at  which  to 

Initiate  run  (sec). 

16-20 

15 

ITIME-Elapsed  number  of  time  steps 
at  which  the  run  Is  Initiated. 

21-25 

15 

ITIMEF-Number  of  time  steps  at 
which  the  run  is  terminated. 

26-30 

15 

ISTART-Restart  control  word: 

0=1nitial  run;  l*restart  using 
previously  stored  flow  solution, 
but  read  Input  from  cards; 

2=restart  using  previously  stored 
flow  solution  and  all  Input  beyond 
card  10. 

3 

1-5 

15 

JAY-Optlon  for  type  of  coordinate 
system:  0=Cartes1an  system  (planar); 
I=cyl1ndr1cal  system  (axl symmetric). 

6-10 

15 

IFLP-Optlon  for  use  of  rotating 
difference  operator:  0=no  rotation 
of  alternating  difference  operator; 
l*rotation  of  difference  operator 
(repeats  every  4  time  steps) 

11-15 

15 

IMESH-Printout  option  flag  for  grid- 
point  output:  0=deletes  printout 
at  all  columns;  l=prints  flow  solu¬ 
tion  at  grid  columns  selected  on 
card  4. 

16-20 

15 

IMESHP-TIme  step  Interval  at  which 
mesh  point  data  will  be  printed 
(when  IMESH=1) 

21-25 

15 

♦INSTREM-Printout  option  flag  for 
flow  solution  along  streamlines: 
0=delete  streamline  output;  l=print 
streamline  output. 

26-30 

15 

*ISTREP-Time  step  interval  at  which 
streamline  output  Is  printed. 

*  It  should  be  noted  that  these  options  have  not  been  fully  coded  for  thermo- 
chemical  output  In  this  version  of  the  program,  and  their  use  Is  discouraged. 


Card 


No. 

Col umns 

Format 

3 

31-35 

15 

36-40 

15 

4 

1-50 

11 

5 

1-5 

15 

6-10 

15 

6 

1-10 

• 

• 

E10. 

• 

71-80 

7 

1-10 

E10. 

71-80 

Description 


♦ISTAT-Printout  option  flag  for 
flow  solution  at  axial  stations: 
0=delete  station  output;  l=print 
station  output. 


*ISTAP-Time  step  interval  at  which 
station  output  is  printed. 


JPRINT-Print  option  flag  array  for 
each  of  the  J  stations  in  the  flow- 
field:  0  placed  in  the  Jth  column 
of  the  card  deletes  the  printout 
for  the  Jth  stations  in  the  grid; 

1  placed  in  the  Jth  column  of  the 
card  allows  for  the  printout  of  the 
Jth  station  in  the  grid. 


NSTRM-Number  of  points  to  be  printed 
along  each  streamline  (maximum=50) 
Note:  delete  card  5  if  ISTREM=0. 


NSL-Number  of  streamlines  to  be 
printed  (maximum=ll).  Note:  delete 
card  5  if  ISTREM=0. 


ZSTRM(I)-Axial  coordinate  of  Ith 
point  along  each  streamline  (1=1, 2, 3. 
...NSTRM);  units  consistent  with  REF 
units.  Note:  delete  card  6  if 
ISTREM=0.  (May  be  continued  on 
additional  cards,  as  necessary.) 


RMF(I)=Ratio  of  weight  flow  rate 
of  Ith  streamtube  to  total  weiqht 
flow  rate  through  the  duct  (1=1,2, 
3....NSL).  e.g.  RMF(1)=0.0,  RMF(2)= 
0.1...RMS(NSL)=1.0.  Note:  delete 
card  7  if  ISTREM-0.  (May  be 
continued  on  additional  cards,  as 
necessary. ) 


] 


*  It  should  be  noted  that  these  options  have  not  been  fully  coded  for  thermo- 
chemical  output  in  this  version  of  the  program  and  their  use  is  discouraged. 
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Card 


No. 

Columns 

Format 

Description 

8 

1-5 

15 

NSTAT-Number  of  axial  stations  at 
which  station  printout  is  desired 
(maximum* 10) .  Note:  delete  card  8 
if  ISTAT=0. 

6-10 

15 

NSTAP-Number  of  radial  points  to  be 
Included  in  the  station  printout 

* 

(maximum=ll).  Note:  delete  card  8 

If  ISTAT=0. 

9 

1-10 

E10.0 

ZSTATH( I)-Axial  coordinates  of  inter¬ 
section  of  axis  or  centerline  and  Ith 

printout  station  ( 1=1, 2, 3. . .NSTAT) ; 
units  consistent  with  REF  units. 

Note:  delete  card  9  if  ISTAT=0. 

(May  be  continued  on  another  card  is 

71-80 

NSTST>8. ) 

10 

1-10 

E10.0 

ZSTATT(I)-Axial  coordinates  of 
intersection  of  wall  point  and  Ith 
printout  station  (1=1, 2, 3... NSTAT); 
units  consistent  with  REF  units. 

Note:  delete  card  10  if  ISTAT=0. 

(May  be  continued  on  another  card  If 
NSTAT>8.) 

71-80 

WHEN  ISTART=2  ALL  THE  FOLLOWING  CARDS  ARE  DELETED 

11 

1-5 

15 

IDOM-Number  of  domains  into  which  the 
duct  is  to  be  divided. 

6-10 

15 

JMAX( I) -Number  of  grid  stations  in  the 
Ith  domain  ( 1=1, 2,3. . . IDOM) 

I  DOM 

Note:  E  JMAX(I)s50-ID0M  is  required. 
1=1 

11-15  15  KMAX( I) -Number  of  radial  grid  points 

in  the  Ith  domain  (I=1,2,3...ID(5m) 
(Maximum=21).  Note  JMAX  and  KMAX 
are  given  in  pairs  for  each  value 
of  I=1,...,ID0M.  e.g.  (JMAX(l), 

KMftX(l)),  ( JMAX(2) ,  KMAX(2) ) .  etc. 


*  It  should  be  noted  that  these  options  have  not  been  fully  coded  for  thermo¬ 
chemical  output  in  this  version  of  the  program  and  their  use  is  discouraged. 
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Card 

No. 

Col umns 

Format 

12 

1-5 

15 

6-10 

15 

11-15 

15 

16-20 

15 

21-25 

15 

26-30 

15 

Description 


INOPT- Inlet  pressure  and  temperature 
option:  0=total  pressure  and  total 
temperature  arrays  at  the  inlet  station 
are  specified;  l=static  pressure  and 
static  temperature  arrays  are  specified; 
lsstatic  pressure  and  static  temper¬ 
ature  arrays  are  specified  at  the  inlet. 


OlttOPT-Discharge  pressure  option: 
0=total  pressure  array  is  specified 
at  discharge  station*, l=static 
pressure  array  is  specified  at 
discharge  station. 

KIN-Number  of  values  to  be  included 
in  the  discharge  station  input  data 
arrays  (maximum=21). 

KOUT-Number  of  values  to  be  included  in 
the  discharge  station  input  data  arrays 
(maximum=21) . 

IRADO-Inlet  radial  velocity  component 
boundary  condition  option :0=specify 
radial  velocity  component,  vr.  (See 
VRINI  on  card  22a);  l=specify  angle 
between  radial  and  axial  velocity 
components,  a  (see  ALFAI  on  card  22b); 
2=specify  ratio  of  circumferential 
component  of  vorticity  to  density, 
n/p  (see  WETAI  on  card  22c). 

ISWIRL-Inlet  circumferential  velocity 
component  boundary  condition  option: 
0=specify  circumferential  velocity 
component,  vs  (see  VTHNIN  on  card 
23a);  l=spec"fy  angle  between  circum¬ 
ferential  and  axial  velocity  compon¬ 
ents.  (See  BETAI  on  card  23b)  '2= 
specify  ratios  of  radial  and  axial 
components  of  vorticity  to  density 
C/p  and  c/p  (see  WXII  and  WZETAI  on 
cards  23c  and  23d).  (This  data  will 
be  ignored  if  JAY=0.) 


*0UT0PT=0  is  not  generally  recommended  since  it  has  been  found  to  become  unstable 
as  the  flow  approaches  a  steady  state. 
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Card 


No. 

Columns 

Format 

Description 

13 

1-5 

15 

NHUB-Nunber  of  coordinate  points 
to  be  specified  along  lower  wall 
axis  or  plane  of  symmetry 
(maximun=30,  minimum=2). 

6-10 

15 

NTIP-Number  of  coordinate  points 
to  be  specified  along  the  upper 
wall  (maximum=30,  minimum=2). 

14 

1-10 

E10.0 

GAM-Ratio  of  spcific  heats,  *y 

11-20 

E10.0 

GASC6N-Gas  constant,  R  ;  dyne 
cm/gnrK 

21-30 

E10.0 

EMINF-Reference  or  plenum  Mach 
number 

31-40 

E10.0 

TTINF-Reference  or  plenum  total 
pressure;  dynes/cm2 

41-50 

E10.0 

PT INF- Reference  or  plenum  total 
pressure;  dynes/cm2 

51-60 

E10.0 

XMDOTF- initial  weight  flow  rate; 
gms/s 

61-70 

E10.0 

REF-Nondimensional izing  refer¬ 
ence  length  (arbitrary  units) 

71-80 

E10.0 

MWINF-Reference  or  plenum  mole¬ 
cular  weiqht;  grams/mole 

15 

1-10 

71-80 

E10.0 

ZBDH(I)-Axial  coordinates  of 
intersections  of  boundary  planes 
with  lower  wall  or  axis  1=1, 2, 3. 
(IDOM+1).  Units  are  consistent 
with  REF.  (For  ID0M=1  ,  1=1 
specifies  beginning  of  domain  1, 
1=2  the  beginning  of  domain  2, 
and  1=3  the  end,  etc.).  May 
be  continued  on  additional  cards 
if  ID(5M>8. 

16 

1-10 

E10.0 

ZBDT(I)-Axial  coordinates  of  in¬ 
tersections  of  boundary  planes 
with  upper  wall  of  duct,  1=1,2.., 
(IDOM+1).  Units  are  consistent 
with  REF.  May  be  continued  on 
additional  cards  if  ID0M>8. 

*Used  only  for  initialization  and  definition  of  entropy  in  output. 
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Card 


No. 

Columns 

Format 

Description 

17 

1-10 

71-80 

E10.0 

RIN(I)-Array  of  radial  coordinates 
of  input  data  points  along  inlet 
station  (1=1 ,2,3. . .KIN) ;  Units  con¬ 
sistent  with  REF.  (May  be  continued 
on  additional  cards  as  necessary.) 

18 

1-10 

71-80 

E10.0 

ROUT ( I ) -Array  of  radial  coordinates 
of  input  data  points  along  discharge 
station  (1-1 ,2,3. .K0UT);  Units  con¬ 
sistent  with  REF.  (May  be  continued 
on  additional,  as  necessary.) 

19a 

1-10 

71-80 

E10.0 

PTINL(I)-Inlet  total  pressure  values 
at  RIN(I)  (1=1 ,2,3. .KIN) , (IN0PT=0) 1 
dynes/cm*.  (May  be  continued  as 
additional  cards,  as  necessary.) 
Cards  19-23  follow  same  format. 

19b 

1-10 

71-80 

E10.0 

PINLETI(I)-Inlet  static  pressure 
array  (1=1 ,2,3. . .KIN) ,  (IN0PT=1); 
dynes/cn2 

20a 

1-10 

• 

71-80 

E10.0 

TTINLI(I)-Inlet  total  temperature 
array  (=2,3. . .KIN) .  (IN0PT=0);  °K 

20b 

1-10 

71-80 

E10.0 

TINLEI(I)-Inlet  static  temperature 

array  (1=1 ,2,3 _ KIN) .  (IN0PT=1); 

°K. 

21 

1-10 

71-80 

E10.0 

VZINI ( I )-Inlet  Axial  velocity 
component  (1=1 ,2,3. . .KIN) ;  cm/sec 
(Only  used  if  local  Mach  number 
is  subsonic.) 

i 


Card 


No. 

Col umns 

Format 

Description 

22a 

1-10 

71-80 

E10.0 

VRINI ( I )-Inl et  radial  velocity 
component  array  (1-1 ,2,3. . .KIN) , 
(IRAD0=O);  cm/sec. 

22b 

1-10 

71-80 

E10.0 

ALFAI(I)-Array  of  inlet  flow 
angle  between  radial  and  axial 
velocity  components  a=tan"l 
(vr/vz),  (1=1, 2, 3. ..KIN 
(IRAD0=1);  radians 

22c 

1-10 

71-80 

E10.0 

WETAI(I)-Array  of  ratio  of 
circumferential  component  of  * 
vorticity  to  density  n/p,  at  the 
inlet  (1=1 ,2, 3.. .KIN),  (IRAD0=2), 
cm^/gm/sec. 

23a 

1-10 

71-80 

E10.0 

VTHINI(I)-Inlet  circumferential 
velocity  component  array  (1=1 ,2,3... 
KIN),  (ISWIRL=0) ;  cm/sec. 

23b 

1-10 

71-80 

E10.0 

BETAI(I)-Array  of  inlet  flow  anqle 
between  circumferential  and  axial 
velocity  components,  B=tan~l  (ve/vz) 
(1=1 ,2, 3.. KIN)  (ISWIRL=1 );  radians 

• 

23c 

1-10 

71-80 

E10.0 

WXII(I)-Array  or  ratio  of  radial 
component  of  vorticity  to  density 
S/p,  at  the  inlet  (1=1 ,2,3. . .KIN) , 

( ISWIRL=2) ;  cm3/gm/sec. 

23d 

1-10 

E10.0 

WZETAI(I)-Array  of  ratio  of  axial 
component  of  vorticity  to  density, 
s/p,  at  the  inlet  (1=1 ,2,3. . .KIN). 
IISWIRL=2);  cm3/gm/sec.  Note: 

WZETAI  data  supplied  when  ISWIRl=2 
only. 

71-80 
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Card 


No. 

Columns 

Format 

Description 

24a 

1-10 

71-80 

E10.0 

PT0UT1 (I ) -Discharge  total  pressure 
values  at  ROUT(I)  (1=1, 2, 3. . .KOUT). 
(0UTT0PT=0);  dynes/cm.  (May  be 
continued  on  additional  cards,  as 
necessary. ) 

24b 

1-10 

71-80 

E10.0 

POUTLI(I)-Discharge  static  pres¬ 
sure  array  (1=1 ,2,3. . .KOUT) . 
(0UT0PT=1);  dynes/cm2 

25 

1-10 

71-80 

E10.0 

TT0UTI(I )-Discnarge  total  tempera¬ 
ture  array  (1=1 ,2,3. .. KOUT)  (only 
required  for  ISWIRL  0)  °K.  (Not 
used.  Zeros  are  permissible.) 

26 

1-10. 

71-80 

E10.0 

ZH(I)-Axial  coordinates  along  the 
lower  wall  axis,  or  plane  of  sym¬ 
metry  (1=1 ,2,3. . . NHUB) .  Units  con¬ 
sistent  with  REF.  (May  be  con¬ 
tinued  on  additional  cards,  as 
necessary. ) 

27 

1-10 

71-80 

E10.0 

RH( I) -Radial  coordinates  corres¬ 
ponding  to  axial  coordinates  ZH 
given  on  card  26  above  along  the 
lower  wall,  axis  (axi symmetric) ,  or 
plane  of  symmetry  (two  dimensional); 
units  consistent  with  REF. 

28 

1-10 

71-80 

E10.0 

ZT ( I ) -Axial  coordinates  of  upper 
wall  (1=1 ,2,3. . . NT IP);  units  con¬ 
sistent  with  REF. 

29 

1-10 

71-80 

E10.0 

RT(I)-Radial  coordinates  of  upper 
wall  corresponding  to  axial  co¬ 
ordinates  Z1  given  on  card  28. 

(1=1 ,2,3. . . NTIP ) ;  units  consistent 
with  REF. 

17 


Columns 


Format 


1-10  E10.0 

71-80 


1-10  E10.0 

■ 

71-80 

1-10  E10.0 

71-80 


_ Description _ 

ALPHAI(i.K) -Array  of  the  Ith  species 
concentrations  in  the  inlet  flow  at 
RIN(K) .  (K=l,2,3. . .KIN) ;  moles/gm 
in  the  current  version  of  R2D2 
the  species  are  input  and  output  in 
the  following  order: 

1=1,  OXYGEN  ATOM  (0) 

1=2,  NITROGEN  ATOM  (N) 

1=3,  ELECTRON  (e~) 

1=4,  ARGON  ATOM  (Ar) 

1=5,  OXYGEN  MOLECULE  (02) 

1=6,  NITROGEN  MOLECULE  (N2) 

1=7,  NITRICOXIDE  MOLECULE  (NO) 

1=8,  NITRICOXIDE  ION  (N0+) 

For  each  species  I  an  array  of  dimen¬ 
sion  KIN  is  read  in  on  a  separate  set 
cards.  i=l,2,3. . .8(max) . 

ALPHAI(9,K)-Array  of  vibrational 
temperatures  for  Oo  in  the  inlet  flow 
(K=l,2,3. . .KIN) ;°K. 

ALPHAI( 10, K) -Array  of  vibrational 
temperatures  for  N,  in  the  inlet  flow 
( K= 1,2,3. . .KIN) ;OKT 

For  other  chemical  systems,  the  first 
S  values  of  the  ALPHAI(1,K)  array 
( 1=1, 2, 3. . .S)  refer  to  the  molar  concen¬ 
trations,  while  the  next  V  values  (v<2S) 
refer  to  the  vibrational  temperatures. 


The  input  for  cards  1  through  29  of  the  equilibrium  version  of  the  program  is  identical 
to  the  nonequilibrium  input.  The  following  input  cards  pertain  only  to  the  equili¬ 
brium  program: 


Card 

No. 

30 


31a 

1-10 

E10.0 

ALPHAI(l.K)- Array  of  first  element 
mass  fractions  in  inlet  flow  at 

# 

RIN(K)  ( K=l,2,3. . .KIN) 

71-80 

31b 

1.10 

E10.0 

ALPHAI(2,K)-Array  of  second  element 
mass  fractions  in  inlet  flow  at 
RIN(K)  (K=1,2,3...KIN).  If  IVG=1, 
this  card  is  omitted. 
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Columns 


Format 


Description 


IVG-Number  of  element  mass  fractions 
to  be  calculated  (maximum=2) 


OUTPUT 


A  sample  of  the  mesh  output  for  the  nonequilibrium  and  the  equilibrium 
versions  is  contained  in  the  Appendix.  Units  are  given  in  the  CGS  system 
corresponding  to  the  assumed  use  of  CGS  units  for  the  input  data. 

The  normal  output  is  broken  into  the  following  categories. 

A  -  Input  data 
B  -  Computed  constants 
C  -  Mesh  point  output 
D  -  Streamline  output 
E  -  Station  output 

Input  Data 

The  input  data  cards  are  printed  in  card-image  format  as  they  are  read,  with 
superimposed  headings  to  identify  each  variable.  The  information  is  therefore 
self-explanatory,  and  is  presented  in  this  manner  to  facilitate  checking  for 
and  correcting  input  errors. 


Computed  Constants 


TINF 

PINF 

RHINF 

AINF 

QINF 

K 

J 

Z 

R 

A,B,C,D 


Reference  or  plenum  static  temperature:  °K 

2 

Reference  or  plenum  static  pressure;  dynes/cm 

■5 

Reference  or  plenum  static  density;  gm/cm 
Reference  or  plenum  sound  speed;  cm/sec 
reference  or  plenum  velocity;  cm/ sec 
grid  row  index 
grid  column  index 

axial  distance;  units  corresponding  to  REF 
radial  distance;  units  corresponding  to  REF 
Components  of  Jacobian  of  transformation  from 
(z,r)  to  (s»n) 


20 


Mesh  Point  Output* 


Z  -  Axial  distance;  units  of  REF 

R  -  Radial  distance;  units  of  REF 

VZ  -  Axial  velocity  component;  cm/sec 

VR  -  Radial  velocity  component;  cm/sec 

VTH  -  Circumferential  velocity  component;  cm/sec 

O 

P  -  Static  pressure;  dynes/crtr 

RHO  -  Static  density;  gm/cm^ 

T  -  Static  temperature  °k 

MACH  -  Mach  number 

MOOT  -  Weight  flow  rate,  j  pv2  +  pv  D/B)irr  dr  evaluated  at 

S  =  constant;  gm/sec 

ENTRO  -  Entropy,  aS/Rq  =  (log(p/pJ  -  Ylog(p/p»))/(Y-l ) 


Block  2  of  the  mesh  point  output  contains  the  species  concentrations  in  moles 
per  unit  mass,  in  the  following  order: 

Oxygen  atom  (0) 

Nitrogen  atom  (N) 

Electron  (e) 

Argon  atom  (Ar) 

Oxygen  molecule  (O2) 

Nitrogen  molecule  (N2) 

Nitricoxide  molecule  (NO) 

Nitricoxide  ion  (N0+) 

Block  3  of  the  mesh  point  output  of  the  nonequilibrium  version  contains  the 
vibrational  temperatures  for  O2  and  N2  in  °K.  For  the  equilibrium  versions 
it  contains  the  element  mass  fractions. 


*This  section  will  be  deleted  if  IMESH=0 


2  2 

Block  4  of  the  mesh  point  output  gives  the  enthalpy  ENTH(cm  /sec  ),  specific 
heats  CPV  and  CVV(cm2/sec2  °K)  .gamma,  GAMF,  and  the  specific  heats  CP  and 
CV(cm2/sec2  °K). 

Streamline  Output  * 

Z  -  Axial  position;  units  of  REF 

R  -  Radial  poistion;  units  of  REF 

V  -  Magnitude  of  the  velocity  vector;  cm/sec 

ALPHA  -  Angle  between  the  radial  and  axial  velocity  components, 

i.e.,  a=tan-1  (vr/vz);  degrees 

2 

P  -  Static  pressure;  dynes/cm 

3 

RHO  -  Static  density;  gm/cm 

ENTRO  -  Entropy,  AS/Rq 

MACH  -  Mach  number 


Station  Output  ** 

The  same  variables  described  above  in  the  Streamline  Output  section  are 
printed  in  the  Station  Output. 


*  This  section  will  be  deleted  if  ISTRM=0 

**  This  section  will  be  deleted  if  ISTAT=0  22 


SYNOPSIS  OF  MAIN  PROGRAM  AND  SUBROUTINES  IN  PROGRAM  R2D2 


A  brief  discussion  of  the  Main  Program,  Principal  Subroutines,  Thermocheml- 
al  Utility  Subroutines  and  the  General  Utility  Subroutines  is  given  for  the 
Program  R2D2.  The  nonequilibrium  and  the  equilibrium  versions  of  the  pro¬ 
gram  differ  only  in  Thermochemical  Utility  Subroutines,  and  both  versions 
are  included  in  the  discussion  thereof. 

The  Thermochemical  Utility  Subroutines  for  the  nonequilibrium  version  and  in 
part  for  the  equilibrium  version  of  R2D2  have  been  adapted  from  a  one-dimen¬ 
sional,  steady,  nonequilibrium  flow  program  developed  by  Calspan  (References 
1,2,3). 


MAIN  PROGRAM 


The  main  program  acts  in  large  part  as  an  executive  routine  which  first  gathers 
together  the  boundary  conditions  and  initial  conditions,  then  advances  the 
solution  through  the  prescribed  number  of  time  steps  by  calling  a  sequence  of 
subroutines,  while  also  calling  a  subroutine  to  print  the  solution  at  pre¬ 
scribed  time  intervals,  and  then  stores  the  final  solution  and  terminates  the 
run. 

A  summary  of  the  major  elements  in  the  logical  construction  of  the  main  program 
are  contained  in  a  flow  chart  presented  as  Figure  (3).  The  principal  compu¬ 
tation  actually  performed  in  the  main  program  is  application  of  the  McCormack 
finite  difference  algorithm  at  the  interior  points.  This  operation  requires 
less  than  20  statements  including  separate  paths  for  ITER=1  and  2  and  separate 
DO  loops  for  the  fluid  mechanical  components  and  the  thermochemical  components 
of  the  system  of  equations.  It  is  identified  in  the  flow  chart  simply  as 
"Compute  VECT's  (predictor  solution)"  and  "Compute  V's  (corrector  solution)." 

A  significant  portion  of  the  coding  preceding  this.operation  is  concerned 
solely  with  definition  of  the  rather  cumbersome  5th  component  of  the  H  array 
associated  with  the  temperature  forms  of  the  energy  equation  (see  Volume  I  of 
this  report). 

Calculation  of  a  series  of  fourth-order  damping  terms  is  included  in  the  coding 
of  the  Main  Program.  However,  the  damping  terms  are  multiplied  by  a  coefficient 
XNU  whose  value  is  define  by  a  DATA  statement  as  0.0.  Since  this  term  has  not 
been  found  to  be  necessary,  it  has  not  been  discussed  elsewhere  in  either  vol¬ 
ume  of  this  report.  However,  the  presence  of  these  terms  is  pointed  out  to 
the  user,  since  they  can  be  activated  by  redefining  the  value  of  XNU  in  the 
DATA  statement.  Values  in  the  range  0.0  <  XNU  <  1.0  may  be  considered  if  a 
need  for  artificial  damping  of  the  solution  appears  warranted.  The  terms  in 
question  appear  only  in  the  momentum  and  energy  equations,  and  are  modeled  after 
the  corresponding  diffusion  and  conductivity  terms. 
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FIGURE  3(a).  FLOW  CHART  OF  PROGRAM  R2D2 

Part  1  of  4:  Initialization. 


•  TIME  -  I  TIME  +  1 
TIME  -  TIME  +  DT 


FIGURE  3(b).  FLOW  CHART  OF  PROGRAM  R2D2 
Part  2  of  4:  Start  of  Time 

Integration  Loop. 
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FIGURE  3(c).  FLOW  CHART  OF  PROGRAM  R202 

Part  3  of  k:  Continuation  of  Time 
Integration  Loop. 
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FLOW  CHART  OF  PROGRAM  R2D2 
Part  k  of  k:  Termination  of  T 
Integration  Loop 


PRINCIPAL  SUBROUTINES 


WALLPOINT  ( ITER.DT) 

ITER  -  iteration  index  for  MacCormack  Algorithm  (ITER=1  for  predictor 
ITER=2  for  corrector  step) 

DT  -  time  step  (Ata^L) 

Performs  the  time  integration  of  the  governing  equations  along  the  duct  walls, 
by  a  modified  version  of  the  interior  point  finite  difference  algorithm,  ex¬ 
cluding  wall  points  at  the  inlet  and  outlet  stations. 

INLET(DT) 

DT  -  tim'-  step  (Ata^/L) 

Imposes  the  Inlet  boundary  conditions  and  performs  the  time  integration  of  the 
equations  at  the  grid  points  spanning  the  inlet  station. 

OUTLET (DT) 

DT  -  time  step  (Ata^/L) 

Imposes  the  outlet  boundary  conditions  (if  the  local  Mach  number  is  subsonic)  and 
performs  the  time  integration  of  the  equations  at  the  grid  points  spanning  the 
outlet  station. 

C0MFGH(R,N,J,K,JAY,IFL0P,J1) 

R  -  radius  array  (r/L) 

N  -  domain  index 

J  -  column  index 

K  -  row  index 

JAY  -  indicator  for  axi symmetric  (JAY=1)  or  planar  flow  (JAY=0) 

IFLOP  -  flag  signaling  end  of  domain  (IFL0P=1  at  interior  columns,  1FL0P=0 

at  first  and  last  columns  of  each  domain). 

J1  -  special  column  index:  Jl=2  for  column  at  which  solution  is  being 

advanced,  Jl=l  for  column  behind  it  and  Jl=3  for  column  ahead  of 

it.  (If  J<3,  J1=J) 
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INDATA 


1.  reads  all  (punched  card)  input  data 

2.  nondimensionalizes  the  variables 

3.  sets  up  finite  difference  grid  (calls  MESH) 

4.  interpolates  inlet  and  outlet  boundary  conditions  arrays  to  grid 
points  on  inlet  and  outlet 

IN  IT 


Initializes  the  flowfield  by  iterating  on  a  value  of  temperature  (or  equivalently 
pressure)  at  each  grid  column  that  yields  a  velocity  which  conserves  mass,  using 
a  constant  total  pressure,  PTINF,  total  temperature,  TTINF,  and  mass  flow  rate 
XMDOTF.  The  angular  velocity  component  initializes  using  conservation  of  angular 
momentum  along  grid  rows,  which  are  reasonable  approximations  to  streamlines. 

The  remaining  components  of  velocity  are  determined  by  linearly  distributing  the 
radial  flow  angle  between  the  angles  at  the  upper  and  lower  walls  (or  axis)  at 
that  grid  column.  Species  concentrations  and  vibrational  energies  are  assumed 
to  remain  frozen  along  grid  rows  (i.e.,  streamlines). 

Note  that  this  initialization  procedure  only  uses  the  following  inlet  boundary 
values:  the  angular  velocity  component,  the  species  concentrations,  and  the 
vibrational  energies.  Thus,  the  remaining  boundary  conditions  are  imposed  im¬ 
pulsively  at  the  first  time  step.  The  values  of  PTINF  and  TTINF  should  be 
selected  judiciously  to  avoid  imposing  an  extreme  impulse  on  the  first  step. 

Subroutine  INIT  also  calls  QIJDOT  to  initialize  the  chemical  and  vibrational  time 
constants  for  subsequent  use  in  determining  the  permissible  time  step. 

OUTPUT (V,V4) 


V  -  three  dimensional  array  of  fluid  mechanical  variables  representing  the 
solution  at  the  current  time  step. 

V4  -  three-dimensional  array  of  chemical  concentrations  and  vibrational 
energies  representing  the  solution  at  the  current  time  step. 
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The  subroutine  prints  the  solution  at  the  specified  time  intervals.  The  de¬ 
pendent  variable  arrays,  V  and  V 4,  are  decoded  into  primative  variables  via  a 
call  to  GETOUT. 

COMFGH 

Computes  the  F,G  and  H  arrays  corresponding  to  the  generic  form  of  the  governing 
equations: 


The  E  array  is  defined  in  the  Main  Program  and  transferred  through  Common  Block 
III .  The  E,F,  G  and  H  arrays  are  triply-dimensioned.  The  first  index  refers  to 
the  variable,  e.g.,  1  for  p ,  2  for  pvz,  etc.  The  second  refers  to  the  column, 
Jl.  The  third  refers  to  the  row,  K. 


COMFGH  is  called  from  the  Main  Program  in  a  loop  which  results  in  definition  of 
E,F,G  and  H  arrays  at  Jl=l,2  and  3  initially,  and  thereafter  only  at  Jl=3,  in 
each  domain.  These  arrays  at  Jl=l  and  2  are  redefined  in  the  Main  Program  after 
completion  of  the  time  integration  for  each  column  by  transferring  the  Jl=2 
column  to  Jl=l  and  Jl=3  to  Jl=2. 


Arrays  E4,  F4,  G4  and  H4  contain  the  chemical  species  and  vibrational  energy 
variables,  whereas  E,F,G  and  H  contain  the  fluid  mechanical  variables. 

AXIS(ITER,  J,N,P,PT) 

ITER  -  iteration  index 
J  -  column  index 

N  -  domain  index 

P  -  pressure  array 

DT  -  time  step  (Ata  /L) 

oo 


Performs  the  time  integration  of  the  governing  equations  along  an  axis  of  sym¬ 
metry. 


STEP(DT) 

DT  -  time  step  (Ata^/L) 

Determines  the  maximum  permissible  size  of  the  time  step  based  on  current  data. 

STREAM(K,  J,  DT,  SDUM,  VRDUM,  RVPDUM,  DVRDT) 

K  -  grid  row  index 

J  -  grid  column  index 

SDUM  -  interpolated  entropy 

VRDUM  -  interpolated  radial  velocity  component 

RVPDUM  -  interpolated  product  of  angular  velocity  component  and  radius 
DVRDT  -  interpolated  time  derivative  of  radial  velocity  component 

Performs  a  stream-path  trace  from  a  grid  point  on  the  outlet  boundary  to  a  point 
on  the  (z,r)  or  (x,y)  plane  at  the  previous  time  step.  Variables  used  in  inte¬ 
gration  of  LaGrangian  form  of  the  equations  are  obtained  by  linear  interpolation. 
May  also  be  called  from  INLET  if  the  axial  velocity  component  becomes  negative 
at  a  point  on  the  inlet  station. 

MESH (I DOM,  ZBDH,  ZBDT) 

IDOM  -  number  of  computational  domains 

ZBDH  -  array  of  locations  of  intersections  of  domain  boundary  stations  with 
the  lower  wall  (or  axis).  Number  of  values  should  be  (IDOM+1). 

ZBDT  -  array  of  Z  locations  of  intersections  of  domain  boundary  stations 
with  upper  wall.  (IDOM+I  values) 

Generates  a  non-orthogonal  grid  which  spans  each  domain  with  a  specified  number 
of  grid  columns  and  rows.  The  first  and  last  columns  are  coincident  with  the 
domain  boundary  stations,  and  interior  columns  are  equally  spaced  between  them. 

The  first  and  last  rows  are  coincident  with  the  lower  wall  (or  axis)  and  the 
upper  wall,  and  the  interior  rows  are  equally  spaced  between  them. 

The  output  from  MESH  is  the  R(J,K)  and  Z(J,K)  arrays  defining  the  grid  point 
coordinates,  the  arrays  of  the  first  derivatives,  DRHDZ  and  DRTDZ,  of  the  lower 
and  upper  walls,  respectively,  and  the  Jacobian  coefficient  matrices,  A,B,C 
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and  D,  which  define  the  transformation  from  (Z,R)  or  (X,Y)  to  computation  co¬ 
ordinates  (c,n). 


GET0UT(J,  K,  DUM,  VR,  VPHE,  VZ,  PI.  RH01,  Tl,  XMACH,  P,  RHO,  ENTRO,  V,  V 4,  ENTH) 


J  -  grid  column  index 

K  -  grid  row  index 

DUM  -  dummy  variable  (not  used) 

VR  -  radial  component  of  velocity  (dimensional) 

VPHE  -  angular 

VZ  -  axial  . 

PI  -  static  pressure  (dimensional) 

RH01  -  static  density  (dimensional) 

Tl  -  static  temperature  (dimensional) 

XMACH  -  Mach  number  (frozen  sound  speed 

P  -  nondimensional  pressure 

RHO  -  "  density 

ENTRO  -  "  entropy 

V  -  array  of  fluid  mechanical  variables 

V4  -  array  of  chemical  concentration  and  vibrational  energy  variables 

ENTH  -  mixture  static  enthalpy  (dimensional) 


Decodes  the  V  and  V4  arrays  into  primitive,  dimensional  variables. 


MASSFL ( DOTM,  I,  VV 1,  VV2,  VV3,  VV5,  VV5) 


2 

DOTM  -  dimensionalizing  factor  for  mass  flow  rate;  L  p  /a 
I  -  column  index 

VV1  -  array  of  interpolated  values  of  p  at  DOTM 

VV 2  -  array  of  interpolated  values  of  pvz  at  DOTM 

VV3  -  array  of  interpolated  values  of  pvr  at  DOTM 

VV 4  -  array  of  interpolated  values  of  pvQr  at  DOTM 

VV5  -  array  of  interpolated  values  of  pT  at  DOTM 


Traces  a  series  of  20  streamlines,  or  lines  of  constant  mass  flow  ration  in  an 
unsteady  flow,  and  interpolates  the  solution  from  the  grid  points  to  points  along 
these  streamlines. 


THERMOCHEMICAL  UTILITY  SUBROUTINES 


SPHEAT  (M,  ALPHA,  CP,  CV) 

W  -  temperature  (T/T  ) 

oo 

ALPHA  -  species  concentrations  (y^MW^) 

CP  -  specific  heat  at  constant  pressure  (Cp/Ro) 

CV  -  specific  heat  at  constant  volume  (Cy/Ro) 

Calculates  specific  heats  (excluding  the  vibrational  energy  partition  function 
for  those  species  specified  as  being  in  vibrational  non-equilibrium).  (Adapted 
from  Reference  1 . ) 

ENTHAL  (W,  A1 ,  A2,  DH,  DE,  ENI) 

W  -  temperature  (T/Tj 

A1  -  species  concentrations  (-y.MW  ) 

I  oo 

A2  -  vibrational  energies  (evi/RoT  ) 

DH  -  mixture  enthalpy  (hMW  /RoT  ) 

DE  -  mixture  internal  energy  (eMW^/RoTj 
ENI  -  species  internal  energies  (e^/ToT^) 


Calculates  the  mixture  enthalpy  and  internal  energy,  h  and  e.  Also  returns  with 
the  species  internal  energies,  ev-.  (The  vibrational  energy  partition  function 
is  excluded  from  the  calculation  for  those  species  specified  as  being  in  vibra¬ 
tional  non-equilibrium.)  (Adapted  from  Reference  1.) 
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QIJDOT  (W,  P,  R,  V,  Al,  A2,  QIJ,  WEVD) 

W  -  temperature  (T/Tj 
P  -  pressure  (p/pj 
R  -  density  (p/p  ) 

CD 

V  -  velocity  (V/a^) 

Al  -  species  concentration  (y.MW  ) 

I  0° 

A2  -  species  vibration  energies  (e^./RoTj 

QIJ  -  time  rate  of  production  of  i^  species  due  to  all  chemical  reactions 
WEVD-  time  rate  of  vibrational  relaxation  of  species 

Calculates  the  chemical  rate  of  production  of  species,  excluding  the  group  from 
i=l  to  e  which  are  calculated  from  conservation  of  elements.  Also  calculates 
the  rate  of  vibrational  relaxation  for  species  in  the  range  from  i  =f+l  to  g, 
which  are  assumed  to  be  in  vibrational  relaxation.  Also  calculates  time  con¬ 
stants  for  chemical  and  vibrational  processes  to  be  used  in  stability  calculation. 
(Adapted  from  Reference  1.) 

WSUM  (IS,  IE,  DELS,  QIJ,  WDSUM) 

I  -  total  number  of  chemical  species 

IE  -  total  number  of  elements  (IE<J) 

DELS  -  coefficient  matrix  for  expressing  WDSUM  in  terms  of  QIJ 

QIJ  -  chemical  rate  of  production  of  species  in  the  range  i=IE+l  to  IS 

WDSUM  -  chemical  rate  of  production  of  species  in  the  range  i =1 , IE 

MOLWT  (ALPHA,  EMM) 

ALPHA  -  species  concentrations  (y^MW^) 

EMM  -  mixture  molecular  weight  (  MW/MW  ) 

oo 

calculates  the  mixture  molecular  weight. 


•. .  - 
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GET  DEL  (IS,  IE,  DELE,  DELS) 

IS  -  nunber  of  chemical  species 
IE  -  number  of  chemical  elements 

DELE  -  two  dimensional  square  array  of  numbers  of  atoms  of  i^ 
element  in  the  species  for  species  1,2, 3,... IE 

DELS  -  two  dimensional  array  of  numbers  of  atoms  of  1^*  element 
in  the  jt(l  species  for  species  IE+1,  IE+2....IS. 

Calculates  the  coefficient  matrix  DELS  used  to  relate  the  rate  of  production 
of  species  i=l  to  IE  to  that  for  species  IE+1  to  IS  using  conservation  of 
elements. 

SUBROUTINE  SUB6  (I,  AIT) 

I  -  species  being  handled  ( 1=1,  species  is  0,  etc.) 

AIT  -  temperature  (T/Tj 

Computes  the  electronic  excitation  contributions  to  enthalpy,  specific  heat 
and  chemical  potential  to  the  Ith  species  (Adapted  from  Reference  1). 
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GENERAL  UTILITY  SUBROUTINES 


MTRXIN  (A,N) 

A  -  matrix  of  dimension  N 
N  -  size  of  matrix 

Inverts  the  A  matrix  and  returns  with  A" *  in  the  same  location. 

IN I 1  (RT,  DUM,  KMAX,  ANS,  R) 

RT  -  radial  position  at  which  value  is  to  be  interpolated 

DUM  -  one-dimensional  array  of  values  of  variable  to  be  interpolated 

KMAX  -  number  of  values  in  DUM  array 

ANS  -  interpolated  value 

R  -  one-dimensional  array  of  radial  coordinates  corresponding  to  the 
DUM  array. 

Performs  a  linear  interpolation  on  a  one-dimensional  array,  and  returns  with 
the  interpolated  value  ANS. 

INT  (RT,  J,  DUM,  KMAX,  ANS,  R) 

Same  as  INTI  except  R  array  is  two-dimensional,  with  J  as  the  value  of  the  first 
index  and  the  second  index  corresponding  to  the  one-dimensional  DUM  array. 

INT2  (RT,  J,  KMAX,  R,  RAT,  IL,  IU) 

RT  -  radial  position  at  which  interpolation  is  to  be  performed 
J  -  value  of  first  index  of  R  array 

KMAX  -  maximum  value  of  the  second  index  of  the  R  array 

R  -  a  two-dimensional  array  of  radial  positions;  the  first  index  is  the 

column  and  the  second  is  the  row. 

RAT  -  ratio  of  distance  from  R(J,IL)  to  RT  to  distance  from  R(J,IL)  to 
R(J,IU) 

IL  -  index  of  row  below  the  position  RT 

IU  -  index  of  row  above  the  position  RT 


Finds  the  indices  of  the  rows  of  the  R  array,  at  column  J,  bounding  the  radial 
position  RT,  and  the  ratio  of  distances  needed  to  perform  an  interpolation. 


SPLINE  (X,  Y,  N,  EM,  SB,  G) 


X  -  array  of  N  coordinates  of  the  set  (X,Y) 

Y  -  array  of  N  coordinates  of  the  set  (X,Y) 

N  -  number  of  points  in  the  X  and  Y  arrays 

2  2 

EM  -  array  of  second  derivatives  of  spline  curve  d  Y/dX  at  each  point 

SB,G  -  arrays  of  parameters  from  spline  curve  (not  needed  for  further  use) 


SPLINT  (X,  Y,  N,  Z ,  MAX,  YINT,  DYDX,  DY2DX,  EM) 


/ 

X  -  array  of  N  coordinates  of  the  set  (X,Y) 

Y  -  array  of  N  coordinates  of  the  set  (X,Y) 

N  -  number  of  points  in  the  X  and  Y  arrays 

Z  -  array  of  X  values  at  which  the  spline  curve  is  to  be  evaluated 

MAX  -  number  of  values  in  the  Z,  YINT,  DYDX  and  DY2DX  arrays 

YINT  -  array  values  of  Y  at  X=Z 

DYDX  -  array  of  values  of  first  derivative  DY/DX  at  X=Z 

DY2DX  -  array  of  values  of  second  derivative  d2Y/dX2  at  X=Z 

EM  -  array  of  second  derivatives  at  the  coordinate  points  (X,Y) 


Subroutines  SPLINE  and  SPLINT  work  in  sequence  as  follows.  SPLINE  fits  a  cubic 
spline  curve  (Y=f(X)  through  the  given  coordinate  points  (X,Y)  and  returns 
with  the  second  derivative  EM=d  Y/dX  at  the  given  coordinate  points,  SPLINT 
may  then  be  called  repeatedly  to  utilize  the  spline  curve  to  interpolate  values 
and  derivatives  of  Y=f ( X )  at  various  points.  The  arrays*  X,Y  in  the  calling 
sequences  of  both  subroutines  should  be  identical,  although  in  principle,  a  sub¬ 
set  of  (X,Y)  can  be  used  in  SPLINT.  The  EM  array  is  output  from  SPLINE  and  in¬ 
put  to  SPLINT,  the  Z  array  is  input  to  SPLINT  and  the  YINT,  DYDX  and  DY2DX 
arrays  are  output.  (Adapted  from  Reference  5.). 


38 


X  -  two-dimensional  array  of  X  coordinates 

Y  -  two-dimensional  array  of  Y  coordinates 

Z  -  two-dimensional  array  of  any  function  of  (X,Y) 

NX  -  initial  guess  of  column  index  (NX>2)  for  search 
NY  -  initial  guess  of  row  index  (NY>2)  for  search 
NDIMX  -  maximum  value  of  first  index  of  X,Y  and  Z  arrays 
NDIMY  -  maximum  value  of  second  index  of  X,Y  and  Z  arrays 
XO  -  X  coordinate  of  point  to  be  interpolated 

YO  -  Y  coordinate  of  point  to  be  interpolated 

ZO  -  linearly  interpolated  value  of  Z  at  X=X0  and  Y=Y0 
I  -  column  index  at  which  search  terminates 
J  -  row  index  at  which  search  terminates 

Performs  a  general,  two-dimensional,  linear  interpolation  of  any  function 
Z-f ( X , Y ) -  The  coordinates  (X,Y)  need  not  be  orthogonal;  any  four  sets  of  points 
should,  however,  form  a  quadrilateral.  The  search  begins  at  indices  (NX, NY)  and 
ends  at  (I,J).  Therefore  repeated  calls  to  LININT  to  interpolate  from  one  co¬ 
ordinate  system  to  another,  for  example,  should  use  NX=I  and  NY=J.  (From 
Reference  5.) 

The  equilibrium  version  of  R2D2,  namely  R2D2EQ,  utilizes  the  same  principal 
subroutines  and  general  utility  subroutines.  The  thermochemical  utility  sub¬ 
routines  have  been  changed  to  calculate  the  equilibrium  chemical  kinetics  as 
prescribed  by  typical  charts  for  equilibrium  air  chemistry. 

The  following  are  Equilibrium  Thermochemical  Subroutines.  The  only  difference 
for  subroutines  SPHEAT,  ENTHAL,  and  SUB6  between  the  equilibrium  and  non¬ 
equilibrium  packages  is  that  the  vibrational  temperature  in  the  equilibrium 
package  is  set  equal  to  the  static  temperature. 


The  following  subroutines  pertain  to  R2D2EQ  only- 

SPHEAT  (VI,  P,  ALPHA,  CP,  CV) 

W  -  temperature  (T/T  ) 

oo 

P  -  pressure  (P/P^) 

ALPHA  -  species  concentrations  (y.MW  ) 

1  00 

CP  -  specific  heat  at  constant  pressure  (C^/ Rq) 

CV  -  specific  heat  at  constant  volume  (Cy/R0) 

Calculates  specific  heats 


ENTHAL(W,  P,  Al,  A2,  DH,  DE,  ENI) 

W  -  temperature  (T/T  ) 

P  -  pressure  (P/P^) 

Al  -  species  concentrations  (y.MHj 

A2  -  dummy  array 

DH  -  mixture  enthalpy  (hMWyR  TJ 

DE  -  dummy  variable 
ENI  -  dummy  array 

Calculates  the  mixture  enthalpy 


P  -  pressure  (P/P^) 

T  -  temperature  (T/TJ 

SUM  -  mixture  molecular  weight  (MW/MWJ 


Calculates  the  equilibrium  mixture  molecular  weight  by  double  interpolation 
of  curve  fits  of  the  molecular  weight  vs  temperature  for  values  of  pressure. 


GAMMA  (Pt  T,  GAM) 

P  -  pressure  (P /Pj) 

T  -  temperature  (T/T  ) 

CD 

GAM  -  equilibrium  value  of  y. 

Calculates  the  equilibrium  y  by  a  double  interpolation  of  curve  fits  of  gamma 
vs  the  temperature  for  values  of  pressure. 


EQUIL  (P,  T,  ALPHA) 

P  -  pressure  (P/P^) 

T  -  temperature  (T/Tj 

ALPHA  -  species  concentration  (y^MUj 

Calculates  the  species  concentrations  by  a  double  interpolation  of  curve  fits 
for  each  species.  The  species  are  curve  fit  vs  the  enthalpy  for  values  of 
pressure.  The  enthalpy  is  obtained  from  a  double  interpolation  of  the  curve 
fits  of  the  enthalpy  vs  temperature  for  values  of  pressure. 
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CHANGES  IN  THERMOCHEMICAL  DATA  AND/OR  PROGRAM  DIMENSIONS 


The  chemical  kinetics  and  thermodynamic  properties  of  the  system  of  reacting 
and  vibrationally  relaxing  gases  are  fully  defined  in  the  BLOCK  DATA  NE  Sub¬ 
routine  of  Program  R2D2NE.  As  mentioned  previously,  the  present  set  of  data 
in  this  subroutine  are  taken  from  the  8  species,  10  reaction  model  of  high 
temperature  air  used  in  References  (1-3).  The  Fortran  names  are  also  largely 
synonymous  with  those  used  in  References  (1-3),  since  the  thermochemical  sub¬ 
routines  employ  portions  of  code  from  the  CALSPAN  Normal  Shock  Program  de¬ 
scribed  therein.  The  block  data  is  defined  in  the  following  dictionary. 


Dictionary  of  Variables  in  BLOCK  DATA  NE 


Name 

Dimension 

Definition 

Present  Vali 

IV 

- 

Number  of  chemical  species 

8 

IVB 

- 

Number  of  elements 

4 

IF,  IG 

- 

Indices  defining  range  of  species 
in  vibrational  nonequilibrium 
(Note:  IF=f+l ,  IG=g) 

5,6 

IVP,  IVBP 

- 

IV+1  and  IVB+1 

9,5 

THETAV(I), 

THEVJP(I) 

8 

Characteristic  vibrational 
temperature  for  i th  species 

See 

listing 

ETAJ(I) 

8 

Number  of  atoms  per  molecule 

4*1.0,  4*2. ( 

SBJ(I) 

8 

Constants  for  chemical  potentials 

See 

listing 

SHJOP(I) 

8 

Heats  of  formation  of  each 
species 

See 

listing 

CNJ(I) 

8 

Number  of  vibrational  levels 
for  each  species 

See 

listing 

CXEF(I) 

8 

Dummy  array  (not  used  -  See 
Reference  1 ) 

See 

1 i sting 

TAUAJ(I) 

TAURJ(I), 

TAUCJ(I), 

TAUDJ(I) 

8 

Constants  which  describe  the 
vibrational  relaxation  rate 
constant  for  species  in  the 
range  f+l<i<q 

See 

listing 

MSUMJ(I) 

8 

Number  of  electronic  levels  for 
each  species  (maximum  of  8  per 
species) 

See 

1 i sting 
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SGJL(I.J) 


8,8 


Degeneracy  of  itfl  electronic 
level  of  jth  species 


CEJLPX(I.J) 

8,8 

Energy  of  ith  electronic  level 
of  j‘n  species 

TFC( I ,J,K) 

5,8,2 

Dummy  array  (not  used) 

CWI(I), 

CZI(I), 

CDI(I) 

10 

Factors,  either  0.0  or  1.0, 
that  control  which  of  3  forms 
of  the  Q. .  relationship  is 

*  J 

used  for  each  reaction 

XNUIJP(I.J) 

3,10 

Stoichiometric  product  coeffi¬ 
cients  for  species  in  ith 
reaction 

XNUI J( I , J) 

8,10 

Stoichiometric  reactant  coeffi¬ 
cients  for  jth  species  in  ith 
reaction 

CAIJ( I ,J) 

8,10 

Indices  identifying  vibration- 
dissociation  coupling  for  jth 
species  in  ith  reaction 

KFIIND(I) 

10 

Indicator  that  denotes  direction 
of  the  ith  reaction  correspond¬ 
ing  to  the  specified  rate  constants 
0  for  forward,  1  for  backword 

AKFI(I) , 
BKFI(I), 
CKFI(I), 

DKFI(I) 

10 

Rate  constants  for  ith  chemical 
reaction:  B.  D. 

k=AiT  1exp(-Ci/T  ') 

FMMLWT ( I ) 

8 

Molecular  weights  of  each  species 

DELE(I,J) 

4,4 

Molecular  weight  of  ith  element 
in  the  jth  species  for  first 

IVB  species 

DELS(I.O) 

4,8 

Molecular  weight  of  ith  element 
in  the  jth  species  for  species 

IVBP  to  IV 

CRO 

- 

Universal  gas  constant 

ISR 

- 

Number  of  reactions 

ISS 

- 

Number  of  species  (same  as  IV) 

See  listing 
See  listing 
80*0 

See  listing 

See  listing 

See  listing 

80*0 

See  listing 

See  listing 

See  listing 
See  listing 

See  listing 

1,98645 

10 

8 


ISF.ISG 

- 

f ,g  (same  as  IF-1 ,IG) 

4,6 

ISC 

- 

Number  of  elements  (same  as  IVB) 

4 

INDSUM 

- 

Indicator  for  inclusion  of  electronic 
excitation  (0-exclude,  1 -include) 

1 

IISFP1 

- 

f+1  (same  as  IF  and  ISFPI) 

5 

IISF,  IISG 

- 

Same  ISF,  ISG 

4,6 

ISCP1 

- 

ISC+1  (same  as  IVBP) 

5 

ISFP1 , ISGP1 

- 

f+1  and  g+1 

5,7 

M 

- 

s+g-f+3  (where  s= I V) 

13 

Ml 

- 

M+l 

14 

MX1 

- 

M-2 

MX2 

- 

M-l 

12 

MX  3 

- 

M 

13 

MX4 

- 

M+l 

14 

IDELXC 

- 

Dummy  variable  (not  used) 

0 

XNUI(I) 

10 

Dummy  array  (not  used) 

10*0 

BBTAI(I) 

10 

Dummy  array  (not  used) 

10*0 

EECHJ(I) 

8 

Dummy  arrays  (not  used) 

3*10*0 

EECNUJ(I), 

EECCPJ(I), 

There  are  obviously  a  number  of  redundancies  in  the  above  list  of  variables, 
which  accrued  during  the  merger  of  the  CALSPAN  routines  with  a  nonreacting 
version  of  the  GASL  duct  flow  program.  These  could  easily  be  eliminated 
by  judicious  use  of  Equivalence  Statements.  Superfluous  data  such  as  M,  Ml, 
MX2,  MX3  and  MX4  could  also  be  calculated,  e.g.,  from  MX1 .  The  dummy  arrays 
have  been  retained  because  they  are,  or  may  have  been,  used  for  temporary 
storage  in  the  thermochemical  subroutines,  or  are  used  under  options  no 
longer  available. 


v/ 


The  dimension  of  8  in  the  thermochemical  data  refers  to  the  number  of  species, 
10  to  the  number  of  reactions,  and  4  to  the  number  of  elements.  These  dimen¬ 
sions  may  be  increased  to  accommodate  larger  systems,  if  necessary,  by 
changing  their  values  in  the  COMMON  blocks  which  appear  in  the  BLOCK  DATA 
Subroutine.  In  addition,  if  the  number  of  species  is  increased  the  dimen¬ 
sions  of  the  dependent  variable  arrays  which  contain  the  species  concentra¬ 
tions  and  vibrational  temperatures,  e.g.,  V4,  VECT4,  ALPHA,  etc.,  must  be 
correspondingly  increased.  Since  every  species  could,  in  principle,  be  a 
diatomic  or  polyatomic  molecule  having  its  own  nonequilibrium  vibrational 
temperature,  these  dependent  variable  arrays  are  dimensioned  to  twice  the 
number  of  species;  16  in  the  present  system.  Virtually  every  subroutine  would 
have  to  be  examined  for  changes  in  the  dimensioned  variables  in  this  case. 

In  some  computations,  it  may  become  necessary  to  increase  the  maximum  grid 
density,  which  is  presently  limited  to  50  columns  and  21  rows.  These  dimen¬ 
sions  appear  in  virtually  all  the  dependent  variable  arrays  as  well  as  the 
grid  coordinate  arrays.  Note  that  the  program  storage  requirement  and  exe¬ 
cution  time  will  escalate  quickly  with  increases  in  grid  density. 

Should  it  be  desirable  to  alter  the  number  of  points  at  which  duct  wall  co¬ 
ordinates  are  specified,  the  dimension  of  30  identifies  the  variables  and 
dimensions  which  will  be  affected.  Similarily,  10  identifies  the  variables 
and  dimensions  associated  with  station  output  and  11  identifies  streamline 
output  variables  and  dimensions.  The  maximum  number  of  points  which  may  be 
printed  along  each  streamline  is  currently  identical  to  the  maximum  number 
of  grid  columns,  50,  and  the  maximum  number  of  points  at  which  inlet  and 
exit  boundary  conditions  can  be  specified  is  identical  to  the  maximum  number 
of  grid  rows,  21.  This  correspondence  in  dimensions  is  coincidental ,  and 
need  not  be  maintained. 

The  dimension  of  5  appearing  in  certain  variables  refers  to  the  number  of 
fluid  mechanical  equations  (continuity,  3  components  of  momentum,  and  energy) 
and  should  therefore  never  be  altered.  The  dimension  of  3  refers  to  the  use 
of  3  grid  columns  in  the  finite-difference  algorithm,  and  is  also  therefore 
invariant.  Similarily  the  dimension  of  2  in  variables  in  Subroutine  WALLPT 
refers  to  the  use  of  2  grid  rows  along  the  walls  of  the  duct,  and  in 
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Subroutines  INLET  and  OUTLET  to  the  use  of  2  grid  columns  along  the  inlet 
and  outlet  planes. 

With  the  above  guidelines  in  mind,  altering  the  program  dimensions  is  a 
straightforward  but  tedious  task. 

Program  R2D2EQ  is  virtually  the  same  as  R2D2NE  except  for  the  replacement 
of  rate  expressions  by  equilibrium  thermochemic^l  properties.  The  data 
contained  in  the  BLOCK  DATA  A  Subroutine  of  R2D2EQ  is  a  subset  of  that 
described  in  the  previous  dictionary.  However,  additional  block  data  is  also 
defined  in  Subroutines  EQUIL,  GAMMA  and  MOLWT,  which  utilize  cubic  spline 
fits  of  tabular  equilibrium  data.  These  subroutines  have  been  tailored,  to 
a  certain  extent,  to  the  considered  air  model.  However,  their  function  is 
clear,  and  they  can  be  easily  recoded  for  a  more  general  model  or  for  other 
specific  models.  For  example,  the  concentrations  of  components  of  the  8  species 
air  model,  the  mixture  molecular  weight  and  the  isentropic  exponent  were  made 
available  (in  graphic  form)  to  the  authors  as  a  function  of  mixture  enthalpy 
and  pressure.  In  addition,  the  mixture  enthalpy  was  available  as  a  function 
of  pressure  and  temperature.  (This  is  not  the  usual  format  in  which  such  in¬ 
formation  is  obtained  but  in  this  case  it  represented  a  readily  available 
equilibrium  solution  which  was  consistent  with  the  nonequilibrium  model.) 
Therefore,  the  procedure  employed  in  these  subroutines  is  to  first  find  the 
enthalpy  at  the  given  pressure  and  temperature,  and  then  find  the  concentra¬ 
tions,  molecular  weight  and  isentropic  exponent  at  that  pressure  and  enthalpy. 
Obviously,  the  enthalpy  is  simply  an  intermediate  variable  in  determining 
these  properties  as  a  function  of  pressure  and  temperature.  Thus,  its  units 
are  irrelevant,  and  it  could  be  eliminated  altogether  if  the  data  were  supplied 
in  a  different  form.  Furthermore,  these  subroutines  use  the  fact  that  one 
of  the  species.  Argon,  is  inert  and  therefore  its  mass  fraction  is  a  constant, 
and  that  two  of  the  species,  NO  and  e",  must  have  identical  mole  fractions 
(to  conserve  charge).  Therefore,  tabular  data  is  only  supplied  for  6  (rather 
than  8)  species. 

The  data  in  Subroutines  EQUIL,  GAWIA  and  MOLWT  is  identified  by  comment  state¬ 
ments  within  each  subroutine.  However,  it  is  reiterated  here  that  there  are 
20  entries  in  the  temperature  table  at  each  of  3  pressures,  or  a  total  of 


60  values.  Thus,  the  temperature  increments  in  the  tables  can  be  different 
at  each  pressure.  There  are  a  corresponding  60  entries  in  each  of  the 
tables  of  enthalpy,  molecular  weight  and  isentropic  exponent,  and  a  corres¬ 
ponding  60  entries  in  each  of  6  species  concentrations  tables.  Since  the 
cubic  spline  fits  can  produce  negative  concentrations  under  certain  con¬ 
ditions,  especially  when  extrapolating,  it  has  been  found  useful  to  fit  the 
logarithm  of  the  concentrations  rather  than  the  value.  An  option  has  been 
provided  for  conversion  of  the  data  to  its  logarithms,  and  to  take  the  antilog 
of  the  fitted  curve,  which  is  controlled  by  ISPEC(I)=0  or  1.  ISPEC(I)  is 
defined  as  1  for  all  species  in  a  DATA  Statement  in  the  present  version  of 
EQUIL,  which  precludes  negative  concentrations. 

Parenthetically,  it  is  also  pointed  out  that  Subroutine  EQUIL  is  only  used 
to  provide  species  concentrations  for  output;  it  is  not  called  by  the  Main 
program  to  carry  out  the  time  integrations  of  the  governing  equations. 
Subroutines  MOLWT  and  GAMMA  are  called  repeatedly  from  the  Main  program  and 
other  subroutines  to  define  the  molecular  weight  and  isentropic  exponent  at 
every  grid  point  on  every  time  step.  In  particular,  the  pressure  must  be 
determined  iteratively  from  the  known  density  and  temperature  using  the  equation 
of  state  and  Subroutine  MOLWT.  It  appears  that  this  process  contributes 
substantially  to  the  execution  of  time  of  R2D2EQ.  Thus  recoding  of  MOLWT  to 
provide  the  molecular  weight  as  a  function  of  density  and  temperature  could 
enhance  the  speed  of  R2D2EQ. 
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APPENDIX 

Sample  Printouts  from  R2D2E0  and  R2D2NE 
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